MOLECULAR CLOUD FORMATION IN SHOCK-COMPRESSED LAYERS 

HIROSHI KOYAMA 2 AND SHU-ICHIRO INUTSUKA 
Division of Theoretical Astrophysics, 
National Astronomical Observatory, Mitaka, Tokyo 181-8588, Japan; 
E-mail addresses: hkoyama@th.nao.ac.jp, inutsuka@th.nao.ac.jp 

0\ ■ ABSTRACT 
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We investigate the propagation of a shock wave into a warm neutral medium and cold neutral 
O [ medium by one-dimensional hydrodynamic calculations with detailed treatment of thermal and 

chemical processes. Our main result shows that thermal instability inside the shock-compressed 
layer produces a geometrically thin, dense layer in which a large amount of hydrogen molecules is 
formed. Linear stability analysis suggests that this thermally-collapsed layer will fragment into 
small molecular cloudlets. We expect that frequent compression due to supernova explosions, 
stellar winds, spiral density waves, etc., in the galaxy make the interstellar medium occupied by 
t> ' these small molecular cloudlets. 
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The tiny scale structures on scales of tens of astronomical units were observed in the interstellar 
medium (ISM). These structures have been detected by 21 cm absorption lines against quasars with VLBI 
techniques (Dieter, Welch, & Romney 1976), against pulsars with time variability (Frail et al. 1994) and 
against close binary stars in optical interstellar lines (Meyer & Blades 1996). The tiny scale structure is 
seen in all directions for which it has been searched, which suggests that it is ubiquitous, not associated 
with large extinction (Heiles 1997). AU-sized structures have been seen not only in atomic gas, but also in 
molecular gas. Langer et al. (1995) observed several clumps ranging in size from 0.007 — 0.021 pc in the 
Taurus molecular cloud 1 (TMC1) with interferometer techniques. The mass of these smallest fragments 
is of the order of O.OIMq. These small-scale structures appear to be gravitationally unbound by a large 
factor. The presence of these tiny scale atomic and molecular structures in variety of the interstellar 
medium suggests that their formation mechanisms have some similarities. 

Models of ISM have been studied by many authors. Field, Goldsmith, & Habing (1969) first showed 
that cold (T ~ 50 K) neutral medium (CNM) and warm (T ~ 8000 K) neutral medium (WNM) can be in 
pressure equilibrium. Cox & Smith (1974) showed that supernova explosions can greatly modify the general 
aspect of the ISM. McKee & Ostriker (1977) developed a theory of a supernova-dominated ISM. The late 
phase supernova remnant (SNR) expands until the internal pressure drops to the ambient interstellar 
pressure (n — 0.1 cm~ 3 , T = 10 4 K). This maximum radius i? m ax = 10 21 pc is reached at a time £ m ax — 
10 6 ' 3 yr. The total volume 5 , i? 3 nax i max which all SNRs occupy in the maximum expansion time t max exceeds 
the volume of the galactic disk, where S ~ 10~ 2 yr^ 1 is the supernova rate in our galaxy. This means that 
the SNRs overlap before they are dissipated. Thus, the study of shock propagation into ISM is important. 

In this paper, we study the basic processes in the evolution of ISM with a strong shock wave by means 
of one-dimensional hydrodynamic calculations with non-equilibrium thermal and chemical processes. The 
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processes we include in our calculation are photoelectric heating from small grains and PAHs, heating and 
ionization by cosmic rays and X-rays, heating by H2 formation and destruction, atomic line cooling from 
hydrogen Lya, C II, O I, Fe II, and Si II, rovibrational line cooling from H 2 and CO, and atomic and 
molecular collisions with grains. Effects of self-gravity and magnetic field are not treated in this paper. In 
the next section, we explain the details of our numerical methods. Results are presented in section 3. In 
section 4, we discuss the formation of small molecular cloudlets. Section 5 is for summary. 



2. NUMERICAL METHODS 

2.1. Hydrodynamic Equations 

To study the thermal and dynamical evolution of shock-compressed layer, we use a one-dimensional 
plane parallel numerical method. The hydrodynamics module of our scheme is based on the second-order 
Godunov method. We use the result of Ricmann problem at each cell interface to calculate numerical 
fluxes, without introducing any artificial viscosity, which characterizes the Godunov method. Higher-order 
Godunov method has been the "state-of-the-art" numerical method in last two decades (see e.g., van 
Leer 1997), and widely used in astrophysical context (see e.g., Truelove et al. 1998) . The hydrodynamic 
equations for the system are written as follows: 

dp dv 
dv DP 

p Tt=-ax> (2) 

P= (l.l+Zo- y)n*BT, (3) 

p = \Amnn, (4) 

1 dP 7 Pdp d 8T 

= ril — nA + — — K— — . (5) 

7-ldt 7-lpdt dX dX y ' 

In the above equations, n is the number density of hydrogen nuclei, p is the mass density of the gas, v 
is the velocity of fluid elements, P is the gas pressure, and K is the coefficient of thermal conductivity. 
We have assumed an abundance of He atom is O.ln. For the range of temperatures and ionized fractions 
considered, the dominant contribution to thermal conductivity is that due to neutral atoms, for which 
K = 2.5 x 10 3 T 1 / 2 cm _1 K _1 s _1 (Parker 1953). We use the ratio of the specific heats 7 = 5/3 for simplicity, 
r and A are heating and cooling rate per hydrogen nucleus, respectively. These rates depend on number 
density n, temperature T, electron fraction x c = n(e~)/n, H2 fraction x% — 2n(H.2}/n, and CO fraction 
^co = fir(CO) /n. The gas behind the shock front is chemical non-equilibrium state. Therefore we solve 
three chemical equations for x e , x%, and xqo and hydrodynamic equations simultaneously. 



2.2. Heating and Cooling Processes 

Wolfirc et al. (1995, hereafter W95) calculated the thermal equilibrium state in the neutral 
atomic phase. We follow their calculation for n < 10 3 cm -3 . We also calculate molecular processes for 
10 3 cm~ 3 < n < 10 6 cm -3 (see Figure [l]). Contributions to the heating of an interstellar cloud are from 
the photoelectric emission from small grains and PAHs, ionization by cosmic rays and soft X-rays, and the 
formation and photodissociation of H2. Following W95, the local FUV field is set to 1.7 times Habing's 
estimate (i.e., Go = 1.7). The cooling function is dominated by line emission from H, C, O, Si, and Fe, by 
rovibrational lines from H2 and CO, and by atomic and molecular collisions with dust grains. To solve 
these thermal processes in non-equilibrium, we solve a set of three time-dependent equations for ionization 
and recombination of hydrogen, and formation and dissociation of H2 and CO. Self-shielding effects are 
also included in the calculation of H 2 photodissociation. The details of these processes are presented in 
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Appendix |A|. Figure [j] shows the equilibrium temperature, pressure, and chemical fractions as functions of 
number density. Important physical timescales are also shown in Figure |l|. 

2.3. Initial and Boundary Conditions 

The strong shock wave is characterized by Mach number. To analyze the evolution of ISM swept up 
by a strong shock wave, we set up a plane-parallel cloud which collides against a rigid wall. Figure |^ shows 
the schematic configuration of our calculation. This configuration also corresponds to a face-on collision 
between two identical clouds. 

We have generated a grid of calculations for varying initial densities and collision velocities. The typical 
models are summarized in table Initial temperatures and chemical fractions are taken to be the value of 
thermal and chemical equilibrium state (see Figure Q). The velocity difference, Vd, defined by upstream 
velocity minus downstream velocity corresponds to the range of Mach number M = 2 — 12. This range of 
Mach number in the WNM (T ~ 8000 K) corresponds to the sweeping-up speed of momentum-driven SNRs. 
External FUV and X-ray radiation penetrate from both sides of two identical clouds. The extinction of the 
radiation is evaluated by the cumulative column density of each grid. The largest total column density we 
considered is that of the standard H I cloud (10 19 ~ 20 cm -2 ). 

3. RESULTS 

3.1. Shock Propagation into WNM 

In this section, we analyze the model W6 where a shock wave propagates into WNM. Figure ^ shows 
the result of calculation. Snapshots at t = 8.0 x 10 4 , 1.3 x 10 5 , 2.5 x 10 5 , and 2.9 x 10 6 yr are shown. The 
horizontal axis in logarithmic scale denotes distance from the rigid wall. Figure [}|a shows the evolution of 
pressure. The shock-compressed layer is almost isobaric. The shock front propagates from left to right. 
Figure |3|b shows the evolution of temperature. Lyman a cooling is so efficient that the temperature of 
post-shock gas quickly decreases to the pre-shock value. In the next subsection, we analyze the shock front 
in detail. Behind the shock front, cooling dominates heating because density is larger and temperature 
decreases monotonically. Around t ~ 10 5 yr, temperature starts to decline rapidly because of thermal 
instability. This instability is mainly driven by O I (63 /im) and C II (157 /im) line cooling. Eventually 
temperature attains thermal equilibrium value T ~ 19 K at n ~ 2.5 x 10 3 cm -3 . Main coolant of collapsed 
layer is C II (see Figure ^jp). Figure shows the evolution of number density. A thin layer is formed by 
thermal instability in shock-compressed layer. The width of thermally-collapsed layer is much shorter than 
the shock-compressed layer. However column density of the thermally-collapsed layer becomes comparable 
after t — 2.6 x 10 5 yr (dotted line on Figure ||c). Figure ||d shows the evolution of velocity. Figure ||e shows 
the evolution of electron fraction. High shock temperature raises the ionization degree in the post-shock 
gas. The timescale of recombination of electron in post-shock layer is larger than the timescale of gas 
cooling. Figure ||f shows the evolution of H2 number fraction. Because the timescale of H2 formation is 
much larger than cooling timescale (see Figure |l]d), this thermal instability is not driven by H2 cooling. In 
the highest density region at final time step (dot-dash), 1.3% of the hydrogen is in H2 and 0.0002% of the 
carbon is in CO. Note that the pre-shock H2 number fraction at final time step (dot-dash) is larger than 
the initial fraction, because H2 photodissociation radiation from the left hand side is shielded by large H2 
column density in the collapsed layer in this one-dimensional plane-parallel model. 
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3.2. Structure Across The Shock Front 

In this subsection, we analyze the detailed structure across shock front which propagates into WNM. 
The same calculation as model W6 at t = 1 X 10 4 are shown in Figure [|. Dashed lines show the calculation 
which used the same initial grid spacing as model W6. We also calculated the model with 320 times higher 
spatial resolution (solid lines) to resolve detailed structure of the shock front. Dotted lines denote the same 
calculation without cooling and heating processes. The width of the shock front is on the order of collision 
mean free path, I = l/na co i ~ 0.0003/rt pc. For a strong shock without cooling, post-shock temperature is 
determined by Rankine-Hugoniot relation as T2 = (7 — l)mVJj /2k&. On the other hands, shock temperature 
with cooling is estimated by the following equation: 

^+E loss = \mVl, (6) 
7 — 1 z 

where £a OS s = / n (A — T) dt w nAL yct (T s )£/Vd. This analytic estimation provides logT s [K] = 4.8 in 
this model. Even the lowest resolution model can provide the same values of physical quantities in the 
post-shock region. 

Radiative precursor which is not treated in our models can change pre-shock temperature and 
ionization. Shull & McKee (1979) have calculated the ionization of the pre-shock gas which is caused by 
the diffuse radiation from the post-shock gas in their models of interstellar shock waves with v s — 40 - 130 
km/s. They have concluded that sufficiently fast shocks (v s > 110 km/s) completely ionize the pre-shock 
gas. To study this effect, we also calculated the model with fully ionized pre-shock gas (model 112), and 
found no significant difference between the model 112 and W12. Therefore we neglect the effect of radiative 
precursor in the rest of this paper. 



3.3. Thermal Instability of Isobarically Cooling WNM 

Figure ^ shows the evolution of the fluid element which has the maximum density at each time 
step. Density, temperature and pressure evolve on the dashed lines from left to right. Solid and thick 
lines correspond to the thermally stable equilibrium. Dotted lines correspond to the thermally unstable 
equilibrium. The pre-shock gas is a thermally stable WNM (n = 0.1 cm~ 3 , T = 8000 K). In Appendix 

, we investigate the stability of isobarically contracting gas by linear perturbation theory to analyze the 
evolution on the dashed lines. We have used the post-shock values of hydrodynamic calculation as the 
chemical compositions in the linear analysis (x e = 0.1, x-i — 10~ 63 , xqq — 10~ 23 ). Because the chemical 
reaction timescale is longer than the cooling timescale. Shaded area in Figure [|b denotes thermally unstable 
region determined by the linear analysis. Larger electron fraction makes cooling rate larger, so that the 
unstable region is wider than the unstable region of equilibrium gas. This shaded region shows that thermal 
instability is inevitable when WNM cools into CNM. 

Next, we discuss about the characteristic length and timescale of the thermal instability. The linear 
growth rate of the thermal instability as a function of perturbation wavelength is calculated with a 
pair of given temperature and density. The growth time of the instability is shown in Figure [| In 
this figure, the vertical axis denotes the unperturbed temperature. The value of the constant pressure, 
Pc/^b = 5 x 10 4 K/cm 3 , is adopted from the result of our non-linear calculation. The density can be 
deduced from the relation, n — P c /kftT. The horizontal axis denotes a quarter of perturbation wavelength. 
Contours of growth time are depicted in this wavelength-temperature plane. In Figure [j]b, the temperature 
evolution of our non-linear calculation (Figure ||b) is superposed upon the growth rate contours (Figure ||a) . 
Thermal instability occurs at temperature less than 7300 K. The instability becomes drastic at temperature 
below 1000 K which corresponds to the growth time less than 10 4 yr. In this way we can understand the 
collapse of the layer in terms of thermal instability. The shortest wavelength of unstable perturbation is 
Amin ~8x 10 _5 pc w 16 AU. The thickness of the thermally-collapsed layer increases with time through the 
accretion of gas. 
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3.4. Shock Propagation into CNM 

In this section, we analyze the model CIO where shock wave propagates into CNM. Although Smith 
(1980) has done similar calculation, his calculation does not have sufficient dynamic range to resolve 
thermally collapsing layer which is of our interest. Figure |?] shows the results of our calculation where 
Mach number M = 10 which corresponds to the velocity difference of 10 km/s. Snapshots at t = 3.1 x 10 3 , 
5.2 x 10 4 , 5.7 x 10 4 , and 8.8 x 10 4 yr are shown. The shock front is almost adiabatic. The post-shock 
pressure can be estimated by Rankinc-Hugoniot relation as follows: 

§ = 1 + 7(7 , + 1 V 2 + Im v /( 7 + 1) 2 M 2 + 16, (7) 

r\ 4 4 

where M is the Mach number in the center of mass frame. When M equals to 10, the post-shock pressure 
is 224 times higher than the ambient pressure. Details of post-shock structure is more complicated than 
that in the WNM case. The post-shock gas, of which temperature and density are determined by adiabatic 
shock condition, is thermally unstable (see Figure ||b). The shock front decelerates after post-shock layer 
collapsed. After thermally collapsed layer is formed, accretion shock is formed inside the shock-compressed 
layer. Accretion shock propagates into compressed layer. In the highest density region at final time step, 
4.9% of the hydrogen is in H 2 and 0.02% of the carbon is in CO. Note that the abundances of these 
molecules are still increasing because the H2 formation timescale is much longer than collapse timescale 
t ~ 10 5 yr (see Figure |Jd). 



3.5. Thermal Instability of Isobarically Cooling CNM 

Figure || shows the evolution of the fluid element which has the maximum density at each time step. 
Shaded area of Figure ||b denotes thermally unstable region. In this case of linear analysis (Appendix 
|b|), chemical compositions in thermal and chemical equilibrium states of each densities are used. This 
Figure shows how thermally stable CNM becomes unstable. The growth time of the instability is shown 
in Figure ^|. In Figure ^|b, we superpose the temperature evolution of our non-linear calculation (Figure 
0b) upon the growth rate contours (Figure ^a). The shortest wavelength of unstable perturbation is 
A min - 2.5 x 10~ 5 pc w 5 AU. 



3.6. Results from Various Initial Parameters 

Calculations are done at nine different initial densities (10~ 2 < n < 10 2 ), and eleven velocity differences 
(2 < Mach number < 12). The column density is 10 20 cm" 2 . We follow the calculation until t = 10 5 - 10 6 



yr. The number fraction of H2 in the highest density region is shown in Figure 10. The horizontal axis 
denotes the initial number density of hydrogen nuclei. The vertical axis denotes the velocity difference. 
Solid lines denote the 8 % boundary of H2 number fraction. Higher velocity produces more abundance 
of H 2 in the thermally-collapsed layer. Note that the fragmentation of the thermally-collapsed layer will 
change the column density of the fragments, which change the efficiency of H 2 self-shielding. Therefore, the 
amount of H2 in the fragments would be different from the value obtained in this paper. To determine the 
realistic abundance of H2, two or three-dimensional hydrodynamic calculation is required. 



4. DISCUSSION 



Our hydrodynamic calculation shows that thermal instability makes the collapsed layer with 
considerable amount of H2. In this section, we discuss the consequence of this instability. 
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4.1. Fragmentation of Thermally-Collapsed Layer 

The shock propagation into ISM can make the thermally-collapsing gas inside the shock-compressed 
layer. The initial thickness of the thermally-collapsed layer is dozens of AU and agrees with the prediction 
by the linear analysis. The thickness of the thermally-collapsed layer increases with time through the 
accretion of gas. The thermally collapsing layer is dynamically unstable also in the Y- and Z-direction, so 
that perturbations in the layer will grow and lead to fragmentation of the layer. We expect that this layer 
will break up into very small cloudlets which have different translational velocities. This velocity dispersion 
of cloudlets produced by the passage of SN shocks may make an origin for the observed interstellar 
"turbulent" velocities. Extremely high pressure in the tiny-scale structure observed in H I clouds (Heiles 
1997) is consistent with our shock propagation model. 

4.2. Molecular Clouds in The Galaxy 

It seems difficult to detect these cloudlets in emission lines of CO molecules (see e.g., Liszt & Wilson 
1993, Liszt 1994). Many shell- like or filamentary structures in the H I 21-cm observation maps (Hartmann 
& Burton 1997) are reminiscent of overlapping shells of old SNRs. If these structures in H I maps are 
really the results of SNR shocks, many tiny molecular cloudlets should be hidden in the shell and should be 
ubiquitous in the Galaxy. McKee & Ostrikcr (1977) showed that the recycling SNRs in the galaxy would 
overlap before they are dissipated. The overlapping filamentary region may produce larger clouds. 

5. SUMMARY 

We have done one-dimensional hydrodynamic calculations for the propagation of a strong shock 
wave into WNM and CNM including detailed thermal and chemical processes. Our results show that the 
shock propagation into WNM can make thin and dense H2 layer by the thermal instability inside the 
shock-compressed layer. The shock propagation into CNM can also make the thermally-collapsed layer. We 
predict that this thermally-collapsed layer will fragment into small molecular cloudlets. Our subsequent 
work which includes two or three-dimensional calculation is aimed to explore how these cloudlets have 
velocity dispersion, and form larger clouds. 

We are grateful to the anonymous referee for helpful comments, which have improved the manuscript. 
We would like to thank Shokcn M. Miyama for discussions and continuous encouragement. 

A. THERMAL PROCESSES 

In this Appendix, thermal and chemical processes in our calculation are summarized. We included the 
formation and cooling of H2 and CO in addition to the processes described by W95. 

A.l. Heating and Cooling Processes 

A. 1.1. Photoelectric Heating from Small Grains and PAHs 

The photoelectric emission from small grains and polycyclic aromatic hydrocarbons (PAHs) induced 
by FUV photons is an important mechanism for heating a diffuse interstellar gas cloud. We use the results 
of Bakes & Tielens (1994, hereafter BT) to calculate the photoelectric heating. In their model, grains are 
distributed with a Mathis, Rumpl, & Nordsieck (1977, hereafter MRN) power-law distribution in size, 
n(a)da cx a~ 3 5 da, between 3 and 100 A in radius. BT calculated the charge distribution for each particle 
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size taking the photoionization and the electron and positive ion recombination into account and found the 
total photoelectric heating from the distribution of particles. The resulting heating rate is given by 

r po = 1.0 x 10~ 24 eG ergs s _1 , (Al) 

where e is the heating efficiency and Go is the incident FUV field normalized to the local interstellar value 
(= 1.6 x 10 -3 ergscm _2 s _1 ) estimated by Habing (1968). BT provide a simple fit to e as a function of 
G TV 2 /n e : 

4.9 x 1CT 2 3.7 x l(T 2 (T/10 4 ) n - 7 

(A2j 



1.0+ [(G T 1 /2/ nc )/i925]0-73 LO+KGoTVa/neJ/SOOO]' 
where n c is the electron density. 

Grains and PAHs may also be an important gas coolant. We treat the heating caused by the electron 



ejection (eq.[ Al]) and the recombination cooling separately, so that the net heating or cooling is the 



difference between the heating and cooling. We use the fit to the recombination cooling provided by BT: 

A po = 4.65 x 10^ 30 r - 94 (G T 1 / 2 /n c )' 3 n o ergs a" 1 , (A3) 

with [3 = 0.74/T - 068 . 



A. 1.2. Ionization Heating by Cosmic Rays and Soft X-Rays 

When the UV flux is strongly attenuated by hydrogen photo-absorption, cosmic rays provide a main 
heating mechanism. We adopt a primary cosmic-ray ionization rate of £cr = 1-8 X 10~ 17 s _1 (including the 
ionization of He) following W95. Primary electrons ejected by cosmic-ray ionization have a larger kinetic 
energy than that of thermal electrons. This primary electron heats surrounding gas by secondary ionization. 
The total heating rate is given by 

r C R = (cRE h (E,x e ), (A4) 

where the function Eh(E,x e ) gives the heat deposited for each primary electron of energy E (Shull & Van 
Steenberg 1985). 

Soft X-rays also ionize and heat ISM. W95 calculated the heating and ionization rates using the 
observed diffuse X-ray spectrum and photo-ionization cross section of H, He, and other trace elements. The 
analytic fits are given in their paper. Soft X-ray is attenuated by photo-absorption of neutral hydrogen 
atom. We use the value of absorbing column density 10 19-20 cm -3 as a standard H I clouds. 



A. 1.3. Formation and Photodissociation Heating from H2 

When a hydrogen molecule is photodissociated the chemical potential of the gas is raised by the 
dissociation energy, 4.48 eV. In addition to this, a small amount of kinetic energy is deposited to the 
photodissociated atoms. According to Hollcnbach & McKee (1979, hereafter HM79), the heating rate is 

r uv = 9i? pump {(2.2eV) [1 + M^/nr 1 } , (A5) 

where 



1.6a;Hexp[-(400/T) 2 ] + 1.4x 2 exp{-[12000/(T+ 1200)1} 



cm" 3 , (A6) 
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and the rate coefficient i? pU mp is given in Appendix A. 2. 2. 

When a hydrogen molecule forms through associative detachment or catalytic reaction, it is in a highly 
excited rovibrational state. Only a small fraction of the formation energy is available for heating the gas, 
while most of the energy escapes the cloud via line emission. We use the heating rates from HM79 as 
follows: 

r gr = i?XH{0.2 + 4.2[l + n cr (H 2 )/7i]- 1 }, (A7) 
and the rate coefficient R is given in Appendix [A.2.2| . 



A. Atomic Line Cooling 

At the high temperatures (T J> 8000 K) collisional excitation of hydrogen Lya can contribute to the 
cooling. We also include low-lying metastable transitions of C II (2326A), O I (6300A), Fe II (5.3 /im), and 
Si II (2240A), using rates from Hollenbach & McKee (1989, hereafter HM89). 

Collisional excitation of the fine-structure lines of C II (158 /im) and O I (63 /im) is a dominant gas 
coolant at temperatures T < 8000 K. Additional cooling is provided by the fine-structure transitions of the 
ground electronic state of Fe II (26 /im), Si II (35 /j,m) (HM89). 

We adopt trace element abundance following W95; xh c = 0.1; xo = 4.6xl0 -4 ; xq — 3xl0 -4 ; xsi = 
3.55xl0~ 6 ; xp c — 7.08xl0 -7 . Several trace elements, C, Si, and Fe, have ionization potentials less than 
13.6 eV and thus are susceptible to ionization by the interstellar UV field. For convenience, the level of 
ionization of these trace elements is assumed to be in singly charged. On the other hands, oxygen remains 
neutral by charge exchange reaction to hydrogen, because oxygen has almost the same ionization potential 
as hydrogen's. 



A. 1.5. H2 Cooling 



At temperatures above 500 K, rovibrational lines of H2 contribute to the cooling function. Following 
HM79 we adopt 

L£(LTE) , _ i'w (LTE) 



A H2 « n(H 2 ) 



1 + n n/n 

where n cr is the critical density defined as follows: 



X2 



(A8) 



,H,Ha 



Ly r ' H2 (LTE) 
L? r ,Ha (n-0)' 



(A9) 



which depends only on temperature, not on density. We use the function 2 (LTE), 2 (n — > 0) from 
HM79. We update the function Lf(n — > 0) from Galli & Palla (1998). We assume an ortho - H2 to para - 
H 2 ratio of 3:1. 



A. 1.6. CO Cooling 

Rotational line emission from CO is a major low-temperature cooling agent for molecular clouds. The 
CO rotational cooling is derived from McKee et al. (1982). The optically thin cooling function is 

, x Mk B T) 2 A , k . 

AC ° TOt = F[U f I TX1 h 1 M/2V A10 

E a [l + {n cr /n) + 1.5(n cr /n)V^j 
where A = 9.7 x 10- 8 s _1 , E /k B = 2.76 K , n CI = 3.3 x 10 6 T 3 3/4 cm- 3 for 12 CO (McKee et al. 1982). 
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The CO vibrational cooling at the lowest level is 

Aco(vib) = 7<n H2 A£io, (All) 

H,H2 

where AEiq — 3080K/fcB- The vibrational rate coefficients for transitions to the v=l state of CO were 
taken from HM89. 



A. 1.7. Atomic and Molecular Collisions with Dust Grains 

The cooling of the gas by cooler dust grains is important at high densities n ;> 10 5 cm' 3 . The cooling 
rate per hydrogen nuclei is 



/ t \ 1/2 /iooA\ 1/2 

A gr = 1.2 x 10- 31 n ij^) ( — ) x[l-0.8exp(-75/r)](r-r gr )ergss- 1 , (A12) 

where T gr is an effective grain temperature, averaged over the assumed MRN size distribution (HM89). 
Since the gas generally collides with the smaller grains, we take T gr to be the temperature 8 K of grains of 
radius a — 100 A. 



A. 2. Chemical Reactions 



For the range of temperatures and densities considered, the recombination time scale is longer than 
the cooling timescale (See Figure |l]i). The residual electrons make the cooling rates larger than that of 
equilibrium abundance. Therefore, it is imparative to treat chemical non-equilibrium processes. 

The main coolant in dense clouds is supposed to be CO, so that CO formation should be taken into 
account. Although H 2 is minor coolant, H 2 enables the formation of CO. Thus, we should take into account 
the formation of H2 . 



A. 2.1. Ionization and Recombination of Hydrogen Atom 



Photoionizing UV fields are attenuated by large absorption cross section of hydrogen atom. The 
hydrogen atom is also ionized by cosmic rays and X-rays. The total ionization rates per hydrogen atom 
are given in the Appendix of W95. Collisional ionization is also important in shock-compressed layers, 
and rate coefficient of this processes is h = 5.8 x 10~ 9 T 4 a50 exp(-15.8/T 4 ) cm- 3 s _1 (HM79). The radiative 
recombination coefficient of hydrogen is given by Shapiro & Kang (1987). 



A. 2.2. Formation and Dissociation of H 2 

H2 is formed catalytically on dust grains and by associative detachment; it is destroyed by collisional 
dissociation, by photodissociation via the Lyman- Werner electronic transitions, and by cosmic rays. 

Following Tielens & Hollenbach (1985, hereafter TH85), the rate at which H2 is produced catalytically 
on dust grains is 

R = 0.5vn H n d o- d S(T) » 6 x 10" 17 (T/300) a5 7i H nS"(T) cnrV 1 , (A13) 

S(T) = [l + 0.04(T + 7 1 gl .)°' 5 + 2 x 10~ 3 T + 8 x 10 _6 T 2 ] _1 , (A14) 

where v is the mean thermal velocity of the atomic hydrogen, nd is the number density of dust grains, and 
ad is the average dust grain cross section. 
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H 2 can be formed by associative detachment of the H ion, 

H + H-^H 2 + e-. (A15) 

We use the formation rates from HM79 to obtain an effective rate coefficient. 

The shock waves can produce sufficiently high temperatures so that collisional dissociation of H2 should 
be considered. The rate of dissociation of H 2 by H atoms is ku — 3.4 x 10~ 9 cxp(8000/T) exp(— 5.19 x 
10 4 /T) cm _3 s _1 (HM79). The rate of dissociation of H 2 by electron impact with LTE vibrational population 
is taken into account (Stibbe & Tennyson 1999). Although these LTE reaction rates are over-estimates for 
H 2 dissociation, there is no significant difference between models with and without collisional dissociation. 
H 2 forms after the gas layer collapses, because the timescale of H 2 formation is much larger than cooling 
timescale. The temperature in collapsed layer is too low for H 2 to dissociate by collision. This is why 
collisional dissociation is ineffective in our models. 

We adopt the rate of photodissociation from TH85, 

i?p Ump = 3.4 x 10- 10 G /?(r) s-\ (A16) 

where 0(t) is the probability that photons of the interstellar radiation field penetrate to optical depth r in 
a plane parallel cloud layer. Photodissociating UV fields are self-shielded by H 2 , and the optical depth r is 



t = 1.2 x 10- 14 A(H 2 )(5T/ D \ (A17) 

where iV(H 2 ) is the column density of H 2 and SVd is the Doppler line width in km s _1 , and we adopt this 
value to be the sound velocity. The analytic expression of 0{t) is given by TH85. 

The rate at which cosmic rays destroy H 2 is estimated by HM89 to be 2.29Ccr, where £cr is the 
primary cosmic-ray ionization rate for atomic hydrogen. 



A. 2. 3. Carbon Chemistry 

We include a simplified treatment of the conversion of singly ionized carbon C + to carbon monoxide 
CO following Langer (1976) and Nelson & Langer (1997). The simplified chemical model assumes the direct 
conversion of C + to CO, or vice versa, without accounting explicitly for the intermediate reactions. The 
equation describing the rate of production and destruction of CO is written 

= k n(C+)np - r co n(CO) (A18) 
where k = 5 x 10 _16 cm 3 s _1 , Too = 10~ 10 Go s _1 , and (5 is defined by the expression 



= KlX[U L > (A19 

P fciz(OI) + G [rcHx/n(H 2 )r 1 

The definition of Tchx is given by Tchx = 5 x 10~ 10 Go, and k\ takes the value k\ = 5 x 10~ 10 . For CO 
photodissociation we do not treat the line self-shielding of CO following Nelson & Langer (1997), because 
we are only interested in the beginning of CO formation. 



B. THERMAL INSTABILITY IN ISOBARICALLY CONTRACTING GAS 

To study the characteristic length and timescale of the thermally-collapsed layer, we have done linear 
stability analysis. Field (1965) analyzed the stability of uniform gas in thermal equilibrium. Schwarz et al. 
(1972) analyzed the isochorically cooling gas without heating. Our non-linear calculation, however, shows 
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that the shock-compressed layer is almost isobaric. In this appendix, we analyze the stability of isobarically 
contracting gas. 



B.l. The Basic Equations 

We consider the gas with number density n, temperature T, velocity v, and pressure P, where all 
variables are functions of position and time. The governing hydrodynamic equations are 

dn dv _ N 

dv dP 

m * n Tt = -ax> (B2) 

1 dP 7 Pdn „ ( d ( T dT\ 

= nr -nA+— [K— , (B3) 



7- 1 dt 7- 1 n dt dX \ dX / 

P = nk B T. (B4) 

In this analysis, we neglect He atoms for simplicity. The effect of chemical evolution can be included in 
non-equilibrium cooling and heating function. For the range of temperatures and densities considered, the 
chemical timescale is longer than the cooling timescale. We use the non-equilibrium chemical fraction just 
after shock heating. 



B.2. The Unperturbed State 

To express isobarically contracting background, we use comoving coordinate x related to real coordinate 

X, 

x = X/a(t), (B5) 

where a(t) is the contraction parameter. In these contracting coordinates, the velocity field can be written 
as 

v = ax + V\, (B6) 
where V\ is the perturbed velocity relative to the contraction. 
Lagrangian differentiation of physical variable / is 



df df ( d\ (df\ if d\ 
d-t^m + { v dx) f= {m) x + a{ V1 ^) f - 

In comoving coordinates, the one-dimensional hydrodynamic equations are 



(B7) 



ct 1 d 

n+-n+- — {nv 1 )=Q, (B8) 
a a ox 



■ a 1 / d \ 1 dP n 

ax + v 1 + -v 1 + - (v 1 —)v 1 + — = 0, B9 

a a \ ox J muna ox 



p+1 a( vi i) p - ^ {" + 1 "} + h - 1)(nA - " r) - {' {9 i 1 
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If we set the perturbed velocity v\ equals to zero in the above equations, the equations governing isobarically 
contracting background are 



a -n = 0, 
a 



.. 1 dPg 

ax H — 

muna ox 

7-+ 7-1 )-7-=r 
a kbJ 



= 0, 



= 0. 



(BH) 
(B12) 

(B13) 



Equation Bll leads n oc a 1 . Temperature is proportional to a because unperturbed gas is isobaric. From 
equation B13, we obtain the approximate solution for a{t) as 



a(t) =a(0)(l -i/7Tnct), 



(B14) 



where r nct is a net cooling time defined by r nc J = (7 — 1)(A - 



T)/k^T. This approximation is valid as far as 

0. 



variations of density and temperature remain small, i.e. < t <C 7T nct . Equation B12 gives dPo/dx 
Thus, we consider uniform unperturbed state, dno/dx = 0. The isobarically cooling gas is characterized 
by a temperature 7b(t) and a number density no(t) which are functions of time, while the pressure Po is 
constant. 



B.3. The Linearized Equations 

We now examine the effect of small perturbations: 



where 



n(x,t) — no(t){X + ni(i) exp(mx)), (B15) 

P(x,t) = P (l + Pi(t)e3tp(iKaf)), (B16) 

T(x,t) = T (t)(l+Ti{t)exp{iKx)), (B17) 

v±(x, t) = v\{t) exp(inx). (B18) 

The first-order equations for T\, Pi, n\, and v\ are 

rii + — v 1 = 0, (B19) 
a 

«! + -vx + c 2 — Pi = 0, (B20) 
a a 

6 . a . , n-y s^Ti + s p ni r T T! + r p rti k 2 

Pi + -Pi - 7«,i H 1 — = -ZCs—Tx, (B21) 

7"net 7"cool 7"hcat ^ 

Pi=?M+Ti, (B22) 

s T = 3 In A/9 In T, (B23) 

s p = 9 In A/9 Inn, (B24) 

r T = d\nT/dlnT, (B25) 

r p = dlnT/dbxn, (B26) 



.-1 



coo 



heat 



l = (7 - l)A/fc B T, (B27) 



(7 - l)r/fc B T. (B28) 
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We adopt the collision mean free path, 1, instead of thermal conduction coefficient K = nk-Qic s /(^f — 1) 
The term a/a remains constant (— l/ 7 T ne t) as long as t <C 7T ne t- We assume that the perturbation is 
proportional to exp(crf) and obtain the dispersion relation as 



.2 2(J , « N2l _ , S T~1 PT-1, ,2 



+ {kc B y > [<r+ — + £c s k z - {kc s ) z I (1 - 7 )<t + — £ 2- ^ = 0, (B29) 

7 r not J \ T coo l Thoat / I T coo i T hoat J 

where fc = n/a. This cubic equation for a has real root and two complex conjugate roots. The former root 
is called condensation mode, and the latter roots are called wave mode. Following Field (1965), we call the 



most unstable mode the condensation mode. In Figure 11, we show an example of dispersion relation of the 
condensation mode. 
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Fig. 1. — (a) Equilibrium temperature and pressure. An absorbing column is 10 19 cm~ 2 (solid lines), 
10 20 cm -2 (dash lines), (b) Chemical fractions as functions of number density, (c) Heating (dashed lines) and 
cooling rates (solid lines) per hydrogen nucleus in equilibrium condition corresponding to panel a (solid lines). 
Heating processes are, photoelectric heating from small grains and PAHs (PE), X-ray (XR), Cosmic-ray (CR), 
and H 2 formation/destruction heating (H 2 ). Cooling processes are CII fine-structure (CII), 01 fine- structure 
(OI), Hydrogen Lyman-a (Ly-a), CO rotation/ vibration line (CO), and Atomic and molecular collisions with 
dust grains (GR). (d) Comparison of various timescales, cooling time (solid line), recombination time (dashed 
line), free fall time (dot-dash line), and H 2 formation time (dotted line) which is defined by (Rnm)^ 1 . 
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Tabic 1. MODELS. 



Model 




T; b 




N H d 


M e 


v d f 


ti s 


n f h 


ly 






W2 


0.1 


8491 


0.03 


1(20)' 


2 


18 


7.0(6) 


1.73 


6353 


3.3(3) 


9.5(-8) 


W6 m 


0.1 


8491 


0.03 


1(20) 


G 


54 


2.9(6) 


2.5(3) 


18.7 


5.2(4) 


0.013 


W6L 


0.1 


8912 


0.05 


1(19) 


6 


54 


2.8(5) 


2.5(3) 


19.0 


5.2(4) 


5.6(-5) 


W12 


0.1 


8491 


0.03 


1(20) 


12 


109 


1.45(6) 


1.2(4) 


16.6 


2.0(5) 


0.21 


112 


0.1 


8491 


1.00 


1(20) 


12 


109 


1.45(6) 


1.2(4) 


16.6 


2.0(5) 


0.21 


C2 


10 


107 


8.7(-4) 


1(20) 


2 


2.1 


7.0(5) 


2.8(2) 


28.5 


8.8(3) 


5.1(-5) 


CIO n 


10 


107 


8.7(-4) 


1(20) 


10 


10 


1.0(5) 


1.4(4) 


16.7 


2.5(5) 


0.049 



a The initial number density of hydrogen nuclei in cm . 
b The initial temperature in K. 
c The initial electron fraction. 

d The total column density of hydrogen nuclei in cm~ 2 . 
e The Mach number in the center of mass frame. 

f The velocity difference in km/s, defined by upstream velocity minus downstream velocity which 
we set to be zero. 

g The final time in year. 

h The maximum density in cm~ 3 at final time step. 

'The temperature in K in the maximum density regions at final time step. 
■'The pressure in K cm~ 3 in the maximum density regions at final time step. 
k The H2 fraction in the maximum density region at final time step. 
'The notation a(b) means a x 10 6 



'Details of this model arc discussed in section 3.1 



"Details of this model are discussed in section 3.4 
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Fig. 2. — Schematic configuration of our calculation. The shock front propagates from left to right. The 
velocity difference, Vd, is defined by upstream velocity minus downstream velocity which we set to be zero. 
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Fig. 3. — The evolution of the shock propagation into the WNM. The figures show (a) pressure, (b) 
temperature, (c) number density of hydrogen nuclei, (d) velocity (e) electron number fraction, and (f) 
H2 number fraction. Solid lines denote the compressed layer at t = 8.0 x 10 4 yr, dashed lines at t — 1.3 x 10 5 
yr, dotted lines at t = 2.5 x 10 5 yr, and dot-dash lines at t = 2.9 x 10 6 yr. 
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Fig. 4. — Structure across the shock front. The model W6 at t = 1 x 10 4 yr is shown. The figures show 
(a) pressure, (b) temperature, (c) number density of hydrogen nuclei, (d) velocity, and (e) electron number 
fraction. Dashed lines correspond to the same spatial resolution as figure ||. Solid lines denote the calculation 
with 320 times higher spatial resolution than dashed lines. Dotted lines denote the calculation with 640 times 
higher resolution without cooling and heating. 
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Fig. 5. — The evolution track of WNM compression on the n - P and n - T plane. Density, temperature and 
pressure evolve on the dashed lines from left to right. Solid and thick lines correspond to the thermally stable 
equilibrium. Dotted lines correspond to the thermally unstable equilibrium. In Appendix |b], we investigate 
the stability of isobarically contracting gas by linear perturbation theory to analyze the evolution on the 
dashed lines, (b) Shaded area denotes thermally unstable region determined by the linear analysis. Larger 
electron fraction makes cooling rate larger, so that the unstable region is wider than the unstable region of 
equilibrium gas. This shaded region predicts that thermal instability is inevitable when WNM cools into 
CNM. 
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Fig. 6. — The growth time of the instability. The vertical axis denotes unperturbed temperature. The value 
of the constant pressure, P c /k B = 5 x 10 4 K/cm 3 , is adopted fr om the result of our non-linear calculation. The 
density can be deduced from the relation, n = P c /k^,T. The horizontal axis denotes a quarter of perturbation 
wavelength. Contours of growth time are depicted in this wavelength-temperature plane. The shortest 
wavelength of unstable perturbation is A m i n ~ 8 x 10~ 5 pc w 16AU. (b) We superpose the temperature 
evolution of our non-linear calculation (Figure ||b) upon the growth rate contours (Figure ||a) . 
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Fig. 7. — The evolution of the shock propagation into the CNM. The figures show (a) pressure, (b) 
temperature, (c) number density of hydrogen nuclei, (d) velocity, (e) electron number fraction, and (f) 
H2 number fraction. Solid lines denote the compressed layer at t = 3.1 x 10 3 yr, dashed lines at t = 5.2 x 10 4 
yr, dotted lines at t = 5.7 x 10 4 yr, and dot-dash lines at t = 8.8 x 10 4 yr. 
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Fig. 8. — The evolution track of CNM compression on the n - P and n - T plane. Density, temperature 
and pressure evolve on the dashed lines from left to right. Solid and thick lines correspond to the thermally 
stable equilibrium. Dotted lines correspond to the thermally unstable equilibrium, (b) Shaded area denotes 
thermally unstable region determined by the linear analysis. This Figure shows how thermally stable CNM 
becomes unstable. 
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Fig. 9. — The growth time of the instability. The vertical axis denotes unperturbed temperature. The value 
of the constant pressure, P c /kB = 2xl0 5 K/cm 3 , is adopted from the result of our non- linear calculation. The 
density can be deduced from the relation, n = Pc/JzbT. The horizontal axis denotes a quarter of perturbation 
wavelength. Contours of growth time are depicted in this wavelength-temperature plane. The shortest 
wavelength of unstable perturbation is A m ; n ~ 2.5 x 10 _5 pc w 5AU. (b) We superpose the temperature 
evolution of our non-linear calculation (Figure |?]b) upon the growth rate contours (Figure ^l) . 




Fig. 10. — The chemical abundance of the thermally-collapsed layers in the one-dimensional hydrodynamic 
calculation. The horizontal axis denotes the initial number density. The vertical axis denotes the velocity 
difference. Solid lines denote the 8 % boundary of H2 number fraction. 
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Fig. 11. — Dispersion relation for the condensation mode of the isobarically contracting uniform gas with 
n = 10 0,25 , T = 10 3 (solid line). The case of thermal equilibrium unperturbed state is also shown (dashed 
line) . 



